library(tidyverse, warn.conflicts = F)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag():    dplyr, stats
library(rvest)
## Loading required package: xml2
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(cluster)
library(ggdendro)
theme_set(theme_light())
source("plota_solucoes_hclust.R")

Tipos de filme de Emma Watson

Emma Watson

Emma Watson

Emma Watson, aka Hermione Granger, eh esta se tornando uma atriz de peso em Hollywood. Diferente de muitos atores mirins ela vem, na fase adulta, alcancando bastante espaco na industria cinematografica. Seus projetos fora da tela do cinema, como o HeForShe, demonstram que ela nao tem medo de se posicionar de acordo com suas opinioes. Sera que essa forca de vontade tambem se reflete nos papeis que ela escolhe retratar?

Nessa analise usarei dados do Rotten Tomatoes sobre os filmes de Emma Watson para encontrar grupos de filmes que tenham caracteristicas semelhantes. Com relacao a renda e nota atribuida pelos usuarios.

from_page <- read_html("https://www.rottentomatoes.com/celebrity/emma_watson/") %>% 
    html_node("#filmographyTbl") %>% 
    html_table(fill=TRUE) %>%
    as.tibble()

filmes = from_page %>% 
    filter(RATING != "No Score Yet", 
           `BOX OFFICE` != "—", 
           CREDIT != "Executive Producer") %>%
    mutate(RATING = as.numeric(gsub("%", "", RATING)), 
           `BOX OFFICE` = as.numeric(gsub("[$|M]", "", `BOX OFFICE`))) %>% 
    filter(`BOX OFFICE` >= 1)

Observando as variaveis separadamente

Uma boa forma de entender os dados com os quais se estaa lidando eh ve-los atraves de alguma visualizacao sumarizada. Os scatterplots abaixo representam os numeros de renda e notas do RottenTomatoes dos filmes.

Rating atribuido no RottenTomatoes

filmes %>% 
    ggplot(aes(x = "Filmes", y = RATING)) + theme_bw() +
    geom_jitter(width = .01, height = 0, size = 2.5, alpha = .6, color="navy") +
    labs(title="Rating dos filmes", x="", y="Rating (%)")

Num primeiro olhar, se pode perceber pelo menos dois grupos formados pelos filmes, um sendo composto exclusivamente pelo filme do rating mais baixo, The Circle de 2017.

Box Office (Renda)

Os valores de box office variam em intervalos maiores que o rating do filme, por isso eh interessante usar alguma tipo de escala para normalizar os dados. Para que valores muito discrepantes nao atrapalhem a definicao dos grupo.

filmes %>% 
    ggplot(aes(x = "Filmes", y = `BOX OFFICE`)) + 
    geom_jitter(width = .02, height = 0, size = 2.5, alpha = .6, color="chartreuse3") +
    labs(title="Box Office dos filmes (sem normalização)", x="", y="Renda (Milhões de dólares)")

A normalizacao utilizada distribui a distancia entre os pontos de acordo com o quanto a distancia representa para a posicao do ponto em questao. Por exemplo, a “distancia” de 1 milhao de dolares significa mais para filmes que tiveram rendimento menor que 5 milhoes do que para filmes que tiveram rendimento maior de 100 milhoes.

filmes %>% 
    ggplot(aes(x = "Filmes", y = `BOX OFFICE`)) + 
    geom_jitter(width = .02, height = 0, size = 2.5, alpha = .6, color="chartreuse4") +
    labs(title="Box Office dos filmes (com normalização)", x="", y="Renda (Milhões de dólares)")

    scale_y_log10()
## <ScaleContinuousPosition>
##  Range:  
##  Limits:    0 --    1

Observando os dois graficos se pode distinguir tres ou mais grupos. O filme The Circle, que teve a menor nota no Rotten Tomatoes, pertence aa massa de pontos localizada no inicio do eixo y, os filmes com os menores box offices. Isso pode ser um indicio de que no conjunto de filmes de Emma Watson os filmes com menores notas no RottenTomatoes tiverem pior rendimento de publico tambem, para validar ou invalidar essa suposicao se utilizam tecnicas de agrupamento.

Agrupamento

Para produzir uma solução de agrupamento precisamos de:

  • Definição de proximidade/distância entre pontos
  • Definição de proximidade/distância entre grupos ou grupos e pontos
  • Processo de agrupamento
  • Decidir quantos grupos existem

Depois vem o principal: avaliar e interpretar a solução. Agrupamento sempre dá um resultado. Nem sempre é útil.

Com uma dimensão

Há duas maneiras principais de agrupar: aglomerativa ou baseada em partição. Vamos explorar primeiro a hierárquica aglomerativa.

distancias.long = filmes %>% 
    select(RATING) %>%
    dist(method = "euclidean") %>% 
    as.matrix %>% 
    reshape2::melt(varnames = c("row", "col"))

distancias.long %>% 
    ggplot(aes(x = row, y = col, fill = value)) + 
    geom_tile()

# distancias = filmes %>% 
#     select(RATING) %>%
#     dist(method = "euclidean") %>% 
#     as.matrix %>% 
#     heatmap()
row.names(filmes) = NULL
agrupamento_h = filmes %>% 
    column_to_rownames("TITLE") %>% # hclust precisa dos rótulos em nomes de linha (ruim)
    select(RATING) %>%
    dist(method = "euclidian") %>% 
    hclust(method = "ward.D")
## Warning: Setting row names on a tibble is deprecated.
ggdendrogram(agrupamento_h, rotate = T, size = 2) 

ggdendrogram(agrupamento_h, rotate = T, size = 2) + 
    geom_hline(yintercept = 45, colour = "red")

Cada junção é um passo do algoritmo. A altura na dendrograma em cada passo significa a dissimilaridade entre os pontos ou grupos juntados naquele passo.

Na medida que vamos aglomerando, as dissimilaridades nas junções tendem a ir aumentando caso haja estrutura de grupos.

data.frame(k = 1:NROW(agrupamento_h$height), 
           height = agrupamento_h$height) %>% 
    ggplot(aes(x = k, y = height)) + 
    geom_line(colour = "grey") + 
    geom_point() + 
    labs(x = "Junções feitas (34 - clusters)", y = "Dissimilaridade na junção")

Vejamos as soluções com diferentes números de grupos.

solucoes = tibble(k = 1:9)

atribuicoes = solucoes %>% 
    group_by(k) %>% 
    do(cbind(filmes, 
             grupo = as.character(cutree(agrupamento_h, .$k)))) 
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
atribuicoes %>% 
    ggplot(aes(x = "Filmes", y = RATING, colour = grupo)) + 
    geom_jitter(width = .02, height = 0, size = 2, alpha = .6) + 
    facet_wrap(~ paste(k, " grupos"))

stats::heatmap() é uma função que visualiza distâncias entre pontos organizando linhas e colunas via hclust:

filmes %>% 
    select(RATING) %>%
    dist(method = "euclidean") %>%
    as.matrix %>%
    heatmap()

Variando o método de linkage

plota_hclusts_1d(filmes, "RATING", 
                 linkage_method = "centroid", # single, complete, average, centroid, median, ...
                 ks = 1:6)
## Warning: Setting row names on a tibble is deprecated.
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

names(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"

Com outras variáveis

Agrupamento sempre dá um resultado. Mesmo quando ele não é útil:

plota_hclusts_1d(filmes, "YEAR", linkage_method = "centroid", ks = 1:6)
## Warning: Setting row names on a tibble is deprecated.
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

Compare as soluções usando a escala linear da variável e a transformada em log:

plota_hclusts_1d(filmes, "`BOX OFFICE`", linkage_method = "centroid", ks = 1:6)
## Warning: Setting row names on a tibble is deprecated.
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

filmes %>% mutate(`BOX OFFICE` = log(`BOX OFFICE`)) %>% 
    plota_hclusts_1d("`BOX OFFICE`", linkage_method = "centroid", ks = 1:6) + 
    scale_y_log10()
## Warning: Setting row names on a tibble is deprecated.

## Warning: Unequal factor levels: coercing to character

Silhouetas

Dada a distância média de um ponto para os demais do seu cluster \(a(i)\) e a distância média do ponto para todos os demais do cluster mais próximo \(b(i)\), a largura da silhoueta de \(i\) é :

\[ s(i) := ( b(i) - a(i) ) / max( a(i), b(i) ) \]

Repare como 1 significa uma boa atribuição para \(i\), 0 significa indefinição e \(-1\) significa que há outro cluster onde \(i\) estaria melhor alocado.

distancias = filmes %>% 
    select(RATING) %>%
    dist(method = "euclidean")

agrupamento_hs = filmes %>% 
    column_to_rownames("TITLE") %>%
    select(RATING) %>%
    dist(method = "euclidean") %>% 
    hclust(method = "complete")
## Warning: Setting row names on a tibble is deprecated.
plot(silhouette(cutree(agrupamento_hs, k = 4), distancias))

plot(silhouette(cutree(agrupamento_hs, k = 2), distancias))

Duas dimensões

p = filmes %>% 
    ggplot(aes(x = RATING, y = `BOX OFFICE`, label = TITLE)) + 
    geom_point() 
p

#ggplotly(p)
agrupamento_h_2d = filmes %>% 
    column_to_rownames("TITLE") %>%
    select(RATING, `BOX OFFICE`) %>%
    dist(method = "euclidean") %>% 
    hclust(method = "centroid")
## Warning: Setting row names on a tibble is deprecated.
ggdendrogram(agrupamento_h_2d, rotate = TRUE)

data.frame(k = NROW(agrupamento_h_2d$height):1, 
           height = agrupamento_h_2d$height) %>% 
    ggplot(aes(x = k, y = height)) + 
    geom_line(colour = "grey") + 
    geom_point() + 
    labs(x = "Número de clusters produzido", y = "Dissimilaridade na junção")

Como sempre, o algoritmo encontra grupos. No caso, parecem até bem separados. Vamos visualizá-los:

plota_hclusts_2d(agrupamento_h_2d, 
                 filmes, 
                 c("RATING", "`BOX OFFICE`"), 
                 linkage_method = "centroid", ks = 1:6)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

O agrupamento está acontecendo todo em função de BOX OFFICE, apenas. Como as escalas são diferentes, BOX OFFICE domina qualquer cálculo de distância euclidiana.

Solução: standardize (aka scale).

agrupamento_h_2d = filmes %>% 
    column_to_rownames("TITLE") %>%
    select(RATING, `BOX OFFICE`) %>% 
    mutate(`BOX OFFICE` = log10(`BOX OFFICE`)) %>% 
    mutate_all(funs(scale)) %>% 
    dist(method = "euclidean") %>% 
    hclust(method = "centroid")
## Warning: Setting row names on a tibble is deprecated.
ggdendrogram(agrupamento_h_2d, rotate = TRUE)

filmes2 = filmes %>% mutate(`BOX OFFICE` = log10(`BOX OFFICE`))
plota_hclusts_2d(agrupamento_h_2d, 
                 filmes2, 
                 c("RATING", "`BOX OFFICE`"), 
                 linkage_method = "ward.D", ks = 1:6) + scale_y_log10()
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

distancias = filmes %>% 
    column_to_rownames("TITLE") %>%
    select(RATING, `BOX OFFICE`) %>% 
    mutate(`BOX OFFICE` = log10(`BOX OFFICE`)) %>% 
    mutate_all(funs(scale)) %>% 
    dist(method = "euclidean")
## Warning: Setting row names on a tibble is deprecated.
plot(silhouette(cutree(agrupamento_h_2d, k = 4), distancias))

Mais variáveis

E se tivéssemos mais de duas variáveis?

filmes2 = agrupamento_h_md = filmes %>% 
    mutate(TITLE_LENGTH = nchar(TITLE)) 

dists = filmes2 %>% 
    column_to_rownames("TITLE") %>%
    mutate(`BOX OFFICE` = log10(`BOX OFFICE`)) %>% 
    select(RATING, `BOX OFFICE`, TITLE_LENGTH, YEAR) %>%
    mutate_all(funs(scale)) %>% 
    dist(method = "euclidean")
## Warning: Setting row names on a tibble is deprecated.
agrupamento_h_md = dists %>% 
    hclust(method = "ward.D")

ggdendrogram(agrupamento_h_md, rotate = T)

cores = RColorBrewer::brewer.pal(4, "Set3")
plot(cluster::silhouette(cutree(agrupamento_h_md, k = 4), dists), col = cores, border = NA)

atribuicoes = tibble(k = 1:5) %>% 
    group_by(k) %>% 
    do(cbind(filmes2, 
             grupo = as.character(cutree(agrupamento_h_md, .$k)))) 
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
atribuicoes_long = atribuicoes %>% 
    mutate(`BOX OFFICE` = scale(log10(`BOX OFFICE`)), 
           YEAR = scale(YEAR), 
           RATING = scale(RATING), 
           TITLE_LENGTH = scale(TITLE_LENGTH)) %>% 
    gather(key = "variavel", value = "valor", -TITLE, -k, -grupo, -CREDIT) 
## Warning: attributes are not identical across measure variables; they will
## be dropped
atribuicoes_long %>% 
    ggplot(aes(x = variavel, y = valor, group = grupo, colour = grupo)) + 
    geom_point(alpha = .4, position = position_dodge(width = .5)) + 
    facet_wrap(~ paste(k, " grupos")) + 
    labs(x = "", y = "z-score")

atribuicoes_long %>% 
    filter(k == 4) %>%
    ggplot(aes(x = variavel, 
               y = valor, 
               colour = grupo)) + 
    geom_boxplot() + 
    geom_point(alpha = .4, position = position_jitter(width = .1)) + 
    facet_wrap(~ grupo) + 
    labs(x = "", y = "z-score")

atribuicoes_long %>% 
    filter(k == 4) %>%
    ggplot(aes(x = variavel, y = valor, group = TITLE, colour = grupo)) + 
    geom_point(alpha = .3, size = .5) + 
    geom_line(alpha = .7) + 
    facet_wrap(~ paste("Grupo ", grupo)) + 
    labs(x = "", y = "z-score")

Analisando a distribuição de box offices dos filmes:

filmes %>% 
    ggplot(aes(x = "Filmes", y = `BOX OFFICE`)) + 
    geom_jitter(width = .01, height = 0, size = 2, alpha = .6)

filmes %>% 
    ggplot(aes(x = `BOX OFFICE`)) + 
    geom_histogram(bins = 16) + 
    geom_rug()

Observando o box office em escala logarítmica:

filmes %>% 
    ggplot(aes(x = "Filmes", y = `BOX OFFICE`)) + 
    geom_jitter(width = .02, height = 0, size = 2, alpha = .6) + 
    scale_y_log10()

filmes %>% 
    ggplot(aes(x = `BOX OFFICE`)) + 
    geom_histogram(bins = 20) + 
    scale_x_log10() + 
    geom_rug()

Relacionando as duas variáveis:

filmes %>% 
    ggplot(aes(x = RATING, y = `BOX OFFICE`, label = TITLE)) + 
    geom_point() 

agrupamento_h_2d = filmes %>% 
    column_to_rownames("TITLE") %>%
    select(RATING, `BOX OFFICE`) %>%
    dist(method = "euclidean") %>% 
    hclust(method = "centroid")
## Warning: Setting row names on a tibble is deprecated.
ggdendrogram(agrupamento_h_2d, rotate = TRUE)

Com escala:

agrupamento_h_2d = filmes %>% 
    column_to_rownames("TITLE") %>%
    select(RATING, `BOX OFFICE`) %>% 
    mutate(`BOX OFFICE` = log10(`BOX OFFICE`)) %>% 
    mutate_all(funs(scale)) %>% 
    dist(method = "euclidean") %>% 
    hclust(method = "centroid")
## Warning: Setting row names on a tibble is deprecated.
ggdendrogram(agrupamento_h_2d, rotate = TRUE)